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ABSTRACT 

The relative gain method of process interaction was 
used for the purpose of developing new control strategies 
for an anaerobic digestion process, using the Andrews and 
Graef’s model. This is done considering both the steady state 
and dynamic approaches. Recent developments like inclusion 
of timelag in growth rate function and inhibition due to 
toxic substances are incorporated in the model. The relative 
gain analysis is done around multiple steady states for 
different hydraulic and toxic loadings, the excess of these 
could result in the failure of a digester. 

The results got from the analysis confirm the 
present day control practice and pave way to explore better 
control possibilities. 



Chapter 1 


INTRODUCTION 


Anaerobic digestion, a traditional biological process 
for treatment and disposal of organic wastes, is receiving 
new and intense attention due to its simplicity and a poten- 
tial to produce combustible gases. In the quest for an alter 
native energy source, a system encompassing energy production 
while controlling the environmental pollution will always be 
preferred. The anaerobic digestion process has gained wide- 
spread acceptance in this respect. 

However, these attributes have been clouded by the 
misconception that anaerobic digestion tends to be a instable 
process. Actually, in many cases of reported instability, 
the difficulties have resulted from faulty operational proce- 
dures. The present day operational practise consists only 
of sets of empirical rules. Hence it was felt that a more 
rational control strategy to put the process operation on a 
quantitative basis, is needed. 

For the purpose of quantifying the process and deve- 
loping a control strategy, the steady state and dynamic rela- 
tive gain method of process interactions is used. The anal- 

4 

ysis is done using the Andrews and Graef's model, where in 
recent developments like inclusion of timelag in growth rate 
function and inhibition due to toxic substances are incorpo- 
rated. The analysis is carried out for different values of 
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the disturbance variables to make the control strategy 
aul table for all types of overloadings, the occurrence of 
which results in the failure of the digester. 
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Chapter 2 

literature survey 

Scientific interest in the gas produced by decomposing 
animal and vegetable materials goes back to the early scien- 
tific revolution in the writings of Robert Boyle and his 
assistant Denis Papin in 1682. Volta in 1776 is given credit 
for first identifying this gas, as he found ’combustible air' 
to be generated everywhere in the neighbourhood of decomposing 
vegetation. William Henry in 1806 showed that this gas was 
identical with the synthetic burning gas, which was later 
called methane. However only in the early twentieth century, 
methanogenesis was found to be connected with microbial acti- 
vity. From then onwards the science of methane production 
during digestion started developing slowly. Today the anaero- 
bic digestion process is considered not only as a possible 
means of recovering energy in the form of methane gas but also 
in reducing the pollution load of organic wastes. The anaero- 
bic digestion process is used widely for treating municipal 
waste sludge and industrial wastes containing high concentra- 
tions of organic materials. It has several significant advan- 
tages over other methods like aerobic digestion, activated 
sludge for the treatment of organic sludges. Other than the 

generation of methane gas, it forms a useful by-product in a 

: 1 
lumus-like slurry that is well suited for land reclamation. 

Jnf 0 r tuna tel y, even with all of these advantages, reports on 
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the application of the process for waste treatment have indi- 
cated difficulties with process stability. 

Actually, in many cases of reported instability, the 
difficulties • have resulted from inadequate or faulty operating 
procedures. ' This is shown by its more successful perfor- 
mance in large waste treatment plants which have technical and 
personnel resources necessary for proper operation. The clas- 
sical approach of the sanitary engineer has been to operate 
the process based on a set of empirical rules and the simplest 
steady state models. As the potentials of the anaerobic diges- 
tion process become more and more evident, there arises a 
great need for a more rational control strategy to put process 
operation on a quantitative basis. An exhaustive evaluation 
of process stability and control strategies for an anaerobic 
digester could entail years of effort if performed on a large 
scale. The tools of system analysis such as modelling and 
simulation, can reduce the time required for such a study and 
can be of considerable value in evaluating the effectiveness 
of different control strategies. They could also be of value 
in improving process design since it would allow comparison 
of the different versions of the process with respect to 
process stability. 

Most models which are used to describe the process 
are steady state models and therefore cannot be used to pre- 
dict process performance during start-up operations or under 
transient conditions resulting from changes in process inputs. 

A dynamic model is called for in this instance, to describe 
the process. 
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Usually the anaerobic digestion process is described 
as a series of reactions, where the insoluble organics are 
converted by the action of extracellular enzymes into soluble 
organics. These are changed by the action of acid producing 
bacteria into volatile acids. These are finally converted 
into methane and carbon dioxide by the methane producing 
bacteria. All these three steps can occur either in one 
reactor or the production of volatile acids can occur in one 
reactor and the methane production in the second reactor. 

2. 1 Development of the Model 

4 

Buswell and coworkers were among the first to study 
the anaerobic decomposition of many organic materials. They 
presented a general formula for the conversion of complex 
organic materials to carbon dioxide and methane. 

C" - I + I V ' tS - I ^ I - 

( 1 ) 

However this formula does not include the fraction of subs- 
trate that is converted to microorganisms, which, though small, 
is necessary for the development of dynamic model of the process, 

A more general formula for the conversion of organic 

5 6 

material to microorganisms, carbondioxide is shown below ' 

Organics ^^^ 4 ^ (2) 

S Cj; Mjj 

where, 

Sv 5= organism yield, moles of organisms produced/mole of 


subs tra te cons umed 
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= carbon dioxide yield, moles of CO^ produced/mole of 
organisms produced. 

Sj^ = methane yield, moles of CH^ produced/mole of organism 
X 

produced. 

The values of and Sjyj can be got from data 

S X X 

of Lawrence and McCarty.^ 

Acetic acid accounts for approximately seventy per- 

7 

cent of the methane produced in anaerobic digestion and the 
dynamic model for the digestion process is for the conversion 
of acetic acid to microorganisms, methane and carbon dioxide. 
Yields will be different for different acids, especially the 
increased ratio of methane to carbon dioxide produced as the 
length of the carbon chain of the acid increases. 


Kinetics of the Process* 

One of the early quantitative expressions for micro- 
bial growth was the Monod equation describing a relationship 
between specific growth rate and substrate concentration 




( 3 ) 


where, 

-1 

JO. = specific growth rate, days 

-1 

~ maximum specific growth rate, day 
Kg = saturation coefficient, moles/liter 
Gg = total substrate concentration, moles/liter. 

However j this Monod function cannot be valid for those subs- 
trates, such as volatile acids, that limit growth at low 
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concentrations and are inhibitory to the organism at higher 

concentrations* Andrews has proposed an inhibition function 

8 

for this purpose, 



(4) 


where, = inhibition coefficient, moles/liter. 

The form in which the substrate exists is also impor- 
tant, and Andrews has further modified tl^e inhibition function 
to consider the unionised volatile acids as the limiting 
substrate.® 

ChS ” + Cj 




1 + 




'HS 


li®. 

^i 


(5) 


where, 

^KS “ substrate concentration, moles/liter 

Gg = ionised substrate concentration, moles/liter 

I 

H = hydrogen ion concentration, moles/liter. 

For pH values above 6, the total substrate concentra- 
tion, Cg, is approximately equal to the ionised substrate 
cohcentration, Gg. Therefore, at a fixed pH and with a known 
total substrate concentration, the unionised substrate concen- 
tration can be calculated from the equilibrium relationship 
for the substrate 
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(H-^) (C") 

^HS ~ ^ 

where, = acid dissociation constant, moles/liter. 

At this stage of development, the model is restricted 
to a constant pH reactor and considers only two variables 
(pH and volatile acid concentration) important for monitoring 
digester operation. This restriction of constant pH was 
removed and the model was extended to incorporate the inter- 
action with bicarbonate alkalinity by considering the carbon 
dioxide-bicarbonate equilibrium. 

(C02)d + H^O ;?== H^ + (HCOg)” 

(H*^) (HCOP ^ ^ 

i. e. ( cn ' — Constant — K. ( 7 ) 

where, 

(C02)jj = dissolved carbon dioxide concentration, moles/liter 

- ionisation constant for bicarbonate, moles/liter. 
The bicarbonate alkalinity in the reactor is related 
to the substrate and net cation concentrations through a 
charge balance as shown in equation (s). For pH > 6, this 
charge balance is simplified to that shown in equations 

(H*^) + (C) (HCOp + 2(003") + (Cg) + (oh") + (A) 

(C) - (A) = (HCOp + (Cg) 

(HCOp = 


or 


( 8 ) 



9 


where, 

C = concentration of cations other than the hydrogen ion, 
eq. /liter 

A = concentration of anions other than those shown in equa>^ 
tion (s), eq. /liter 

C 2 = cation concentration, (G - A), eq. /liter. 

With the bicarbonate concentration now known from 
equation (s), the pH is calculated from equation (?) for cons- 
tant values of dissolved carbon dioxide. 

Equation (?) can be rewritten as 


(H"^) (C2 - Cg) 

5^ = ''1 


(9) 


where, = (CO^)^ = dissolved carbon dioxide concentration, 
moles/liter. 

From equations (6) and (9), we get a relation for 
unionised substrate concentration 


-Cg 

'HS ~ 


(10) 


The restriction of constant dissolved carbon dioxide 
concentration was removed by considering the interaction bet- 
ween the gas and liquid phases of the reactor. At equilibrium, 
the dissolved carbon dioxide concentration is related to the 
partial pressure of carbon dioxide in the gas phase by Henry’s 
law 


(OOg)^ 


(co^)^ 
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(CO^)^ = KP^ 


( 11 ) 


concentration of carbon dioxide in gas phase, 
moles /liter of gas volume 

concentration of dissolved carbon dioxide in the 
liquid phase when in equilibrium with the gas 
phase, moles /I iter 

partial pressure of carbon dioxide in gas phase, atm. 
Henry’s law constant, 0.024593 moles/atm-liter 
at 38"C. 

Under normal operating conditions, the carbon dioxide 
in the gas and liquid phases can be considered to be in equi- 
librium. However, this was not so under dynamic conditions 
since transfer across an interface is slow compared with 
ionic reaction rates. An expression for gas transfer was 
incorporated into the model to increase its applicability. 

This expression calculates the rate of transfer of carbon 
dioxide either to and from the liquid phase. 

= la CKPc 

where, 

Tg = gas transfer rate, moles /liter-day 

-1 

gas transfer coefficient, day . 

Andrews and Graef found that in addition tof hydraulic 
and organic loading, process failure can also be caused by 



where, 

(C02)q •= 

Pc 

K 
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1 0 

toxic materials in the feed stream to a digester. As a 
first approximation they assumed that the rate of organism 
kill is first order with respect to the concentration of the 
toxic agent as defined in equation 


■K 


K.J. 


(13) 


where, 

= rate of organism kill, moles/liter/day 
Kj = toxicity rate coefficient, moles of organism killed/ 
mole of toxic material/da y 
Cj = concentration of toxic material, moles/liter. 

Kubicek et al. studied the effects of timelag in the 
specific growth function and incorporated this effect as 
the change in specific growth rate in the system of model 
equation. This was done as the microorganisms present in the 
sludge do not react immediately on the change of the conditions 
in the environment. The increase of inhibition by the subs- 
trate can also be described as increase of the timelag (delay) 
in the growth function. Thus we have 




„ /r. ^ 


K. 

1 


(14) 


Here Gj^(t ‘"1^) denotes the concentration of the unionised 
substrate at the time (t - tj^) and time t^^ describes the 
magnitude of the timelag. 


2.2 Stability of the Digester 

Instability that leads to the ultimate failure of a 
digester can result from any material or effect that interferes 
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with the methane formation. The potential causes of digester 
instability are hydraulic, organic and toxic overloading. 

This is discussed in detail in Chapter 3. 

2.2.1 P rocess Condition Indicators 

Simulation studies have evaluated the potential 
process condition indicators using this three types of over- 
loading. 

Total Volatile Acids ? 

The total volatile acids (TVA) concentration has been 

monitored in most digesters as an indication of the process 

condition of the digester. The simulations and reports from 
1 1 

literature indicate that a sharp rise in TVA can serve as a 

warning and suggest that one TVA analysis per digester per 

day should be adequate for routine surveillance. Nowadays, 

12 

gas-liquid chromatography thwarts the difficulties in 
on-line measurement of TVA concentration and can expedite 
this TVA analysis making TVA a more useful process condition 
indicator. 

Rate of Methane Production* 

This is a process variable which has not been widely 

13 

used. The theory of the digestion process indicates that 
process condition of the digester is directly related to the 
rate of methane production* In all the three cases of loadings, 
the rate drops sharply as failure ensues, particularly in the 
case of toxic overloading where the rate starts to decline 
immediately. It is this characteristic which might enable 
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to sense immediately failure due to a foreign toxic subs- 
tance. Another useful feature of this variable is that it 
can be calculated from readily available on-line gas phase 
measurements. 

Carbon Dioxide^ 

1 4 

Guarino et al, has reported the interest in the 
carbon dioxide content of the digester off gas for the purpose 
of monitoring the process. Their simulation study indicates 
that the carbon dioxide, composition can fluctuate when the 
three types of loadings approach failure conditions. Carbon 
dioxide, the end product of the digestion process partially 
escapes into the gas phase and the tjalance either remains in 
solution as dissolved carbon dioxide or reacts with a base 
such as ammonia to form bicarbonate ions. The solubility of 
carbon dioxide is a function of partial pressure and the 
conversion of carbon dioxide to carbonates and bicarbonates 
depends upon partial pressure and pH. 

pHJ 

Literature indicates that pH is a process variable 
that has been monitored in many digester installations to 
supplement infoimation got from other process variables. 
Simulation studies indicate that the pH declines as the 
process begins to fail. This decline results from the accu- 
mulation of fatty acids and should |>rovide an indication of 
changes in the condition of the process. 
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Other Process Variables s 

Other variables that have been monitored in digestion 
practise to warn of process failure include total alkalinity, 
bicarbonate alkalinity and rate of total gas production. 

None of these variables were effective in predicting digester 
failure. However bicarbonate alkalinity would be a good 
indicator of digestion condition due to its relationship with 

the volatile acids concentration. 

1 5 

Kubicek et al. indicated that the amount of micro- 
organisms would be a good indicator of the digester condition 
and this can either be measured or the microorganism activity 
can be estimated from the methane gas production rate. 

2. 2 Operational Factors Affecting Stability 

1 6 

Andrews and Graef identified several factors that 
would add to digester stability against overloading. These 
included residence time, organism recycle, bicarbonate alka- 
linity, influent substrate concentration and frequency of 
loading. 

Residence Time* 

Increasing the residence time improves digester 
stability against overloading because the incoming substrate 
is immediately diluted to a lower concentration and provides 
the organisms with additional time to metabolize the subs- 
trate. The effective residence time can be increased by 
increasing the reactor volume and/or decreasing the reactor 
flow rate. The effective volume can be increased by having 
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additional digesters, by improving mixing to remove dead 
zones in the reactor or by removing grit and scum accumula- 
tions in the digester. 

Bicarbonate Alkalinity! 

16 

Simulation results suggest that increasing bicar- 
bonate alkalinity enables the digester to sustain sizable 
step increases in step loading. One way to increase bicar- 
bonate alkalinity is by adding a base. Another way, is by 
increasing the influent sludge concentration to the digester. 
However, a sludge that is too highly concentrated can, 

through protein degradation, release inhibitory levels of 

17 18 

ammonia. * So most sludges should be concentrated to 
increase bicarbonate alkalinity without yielding detrimental 
concentrations of ammonia. 

Influent Substrate Concentration! 

Simulation results indicated that a concentrated 
influent sludge is more beneficial to digester stability than 
a dilute sludge. A concentrated sludge feed yields a higher 
steady state organism concentration within the digester, 
thus enabling it to sustain a larger increase in organic 
loading* 

Digested Sludge Return! 

Stability can be gained by providing sludge recycle 

for digesteis. This is usually done by various methods — 

19 

recycling a portion of digested sludge to the reactor, by 
inoculating a digester with well digested sludge from another 
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reactor or by degasifying and concentrating sludge in a 

20 

solid-liquid separator followed by return to the digester. 
Sludge return serves to inoculate the digester with additi- 
onal viable organisms. Thus when overload occurs the added 
organism serve to increase the rate of substrate decomposition. 

Loading Frequency^ 

1 6 

Studies . indicate that stability decreases as the 
frequency of loadings decrease. Certain digesters are loaded 
on a batch feed or discontinuous basis. The larger oscilla- 
tions of the simulation variables can be substantially 
reduced by loading the digesters several times per day as 
compared to loading them once a day. 

2. 3 Control Strategies 

Much of the earlier work on digester control has been 
directed towards preventing or off-setting failure. Little 
work has been done toward using process control techniques 
to increase the organic loading rate that a digester can 
handle oc to improve the final quality of digester effluent. 
Control actions available in anaerobic digestion include a 
decrease in the rate of organic loading, a temporary halt 
in organic loading, addition of a base such as lime and 
soda a&h, dilution of the digester contents or the addition 
of well digested sludge from another digester. Andrews and 
Graef suggested a control action which includes scrubbing of 
carbon dioxide from the digester gas and subsequently recy- 
cling the gas. This provides process control by adjusting 
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the pH of the digester through the removal of a weak carbonic 
acid, in contrast to the usual practise of addition of a base. 
pH is the measured variable and when it deviates from its 
set point, all on-off controller diverts a portion of the 
gas into a gas scrubber where carbon dioxide is eliminated by 
absorption in a- water spray. The scrubbed gas is then recir- 
culated to the reactor, thus reducing the partial pressure 
of carbon dioxide in the gas phase. This, in turn, causes 
dissolved carbon dioxide to be vented from the solution 
phase, giving an increase in pH due to the removal of 
carbonic acid. 

Another control strategy for prevention of failure 

by overloading of toxic materials is the recycling of the 

concentrated sludge from a second stage digester using the 

1 3 

rate of methane production as a feedback signal. The rate 
of methane production can be easily calculated and would be 
an excellent indication of the activity of the methane bac- 
teria that are the most sensitive and critical organisms in 
the digester. 

These current control strategies of the digester 
operation are based on empirical rules evolved from practical 
experience. Kubicek^^ suggested a need for a control system 
which will circumvent indirect regulation methods without 
being too sophisticated to remain practical. 
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Chapter 3 

CHARACTERISTICS OF ANAEROBIC DIGESTION PROCESS 


3 . 1 Dynamic ^'^del of the Process 

The final equations representing the dynamic model 

for the process are now presented by approximating anaerobic 

digesters as continuous flow, stirred tanl< reactors (CFSTR) 

and by making material balances on the appropriate components. 

The general form of a material balance is as follows i 

Rate of material Rate of appearance or 

flow into reactor disappearance of material 

due to reaction 


Rate of material flow Rate of accumulation of 
out of reactor material in reactor 


For a multi -phase process such as anaerobic digestion 
it may be necessary to make material balances for specific 
components in each phase as well as to consider transfer of 
the components between phases. 


(a) The mass balance for the organism 

dCjj 


V 


= F(Cx^ - Cj,) + - K^C^.V 


(15) 


where, 

F = liquid flow rate, liters/day 

V = liquid volume of reactor, liters 

0 — subscript indicating influent concentration. 

(b) The material balance for the total substrate is* 
dCc ' 



20 iq 


(c) Mass balance on carbon dioxide dissolved in the solu- 
tion phase 

dCp 

V = F(G^ - C^) + R^.V + R^.V + T^.V 


dt 


'Cq -C' ‘ “B** ■ “C" ■ *G- 




‘^X ^C^ 


(17) 


F r , dCe dC, 

V + dT ■ dT 


(d) Gas phase carbon dioxide balance 


dPp 

'dt' 


= -PS 


)L. T C Q 

P Vq * ^G Vq 


(18) 


(e) Net cation balance 


dC- 

'"dFf = 


(f) Net toxic chemicals balance 
dC_ 


(19) 


( 20 ) 


The model of a digester and the information flow are 
indicated in Figure 3. 1. • 


3»2 Parametric Uncertainties 

The dynamic model of the digestion process consists 
of several parameters like and K. 

There exists an uncertainty in the value of certain para- 
meter like and because they are difficult to be 

estimated. Other parameters are ’certain’ or their values 
are known precisely. Uncertainty arises due to the conditions 
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of the environment-temperature, concentration and composi- 
tion of substrate. 

Saturation Constant 

There is an uncertainty in the value of Kg 
to be expressed in concentrations of volatile fatty acids 
but it is now expressed in terms of the limiting substrate, 
the unionised acid. This unionised acid is dependent on an 
accurate knowledge of pH. Therefore Kg is very difficult to 
be estim.ated. The value of saturation constant used is 

Kg = 0.333 X 10“^ mol/lit at 38*C.^^ 

Gas Transfer Coefficient, 

The values of gas transfer coefficient are not avail- 
able for anaerobic digester. The model representing the 
system is stiff for high values of However, for the 

values of > "lO* "the results practically do not depend 

on the value of Kj^^^ as the liquid and the gas phases are in 

1 5 

equilibrium. The value of used here is 30/day. 

Inhibition Constant, 

Like there is an uncertainty in the value of 
K^. K^ used to be expressed in concentrations of volatile 
fatty acids. Now, unionized fatty acids are found to be 
the limiting substrate and hence is expressed in terms of 
unionised fatty acids. . The unionised fatty acids are 
very difficult to be estimated and this causes the uncer- 
tainty in K^. is dependent on temperature, concentration 


K« used 
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and composition of the substrate. This dependence, however, 

is not yet quantified. Although values of for different 

volatile acids are not available, there is some evidence that 

some volatile acids are more inhibitory than others. Buswell 
2.1 

and Morgan showed that a digester works better with a 
single substrate acetic acid than with both acetic and pro- 
pionic acid. Venkobachar indicates that by properly moni- 
toring the conditions of the environment a higher fatty acid 
like lactic acid can not only be less inhibitory but can also 
lead to higher methane production rates. 

The value of IC used here is 1.0 x 10~^ mol/lit. 


Ionisation Constant for Fatty Acids, 

The value of K is certain and can be determined, 
a 

depends upon temperature, coinposition and concentration 

of the substrate. However, for lower fatty acids does 

not seem to change appreciably with temperature, were Xo. 
o adCo X '6"^ r»\c*es/ab . 

Ionisation Constant for Bicarbonate, * 

The value of is certain and can be determined. 
The value of used here is 0.1 x 10 mol/lit. 


Henry’s Law Constant, 

Henry's law constant K, too is certain and can be 
determined. The value of ^ used here is 0.024593 mol/atm-lit. 

The average value of the parameters Kg, K^, and 
yields Sj^ and for acetic, propionic and butyric acid 
are given in Table 3.1* 

. L 
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Table 3,1 


Average Values of 
and Butyric Acids 

the Parameters 

for Acetic, 

Propionic 


Parameter 

Acetic acid 

Propionic acid Butyric ac 

(1) 

Ke (mg/lit) 

* 35 “C 

154 

32 

5 

(2) 

K. (mol/lit) 

^ 38 “C 

1.0 X 10"^ 

- 

- 

(3) 

K. (mol/lit) 

35‘‘C 

^0-4.56 

lo-*^.^ 

1 

(4) 


47.0 

47.0 

80.0 

(5) 


47.0 

47.0 

47.0 


The values of 

parameters used by Andrews 

and Graef 


are given in Table 3.2i 

3. 3 Multiplicity of Solutions 

Andrews and Graef in their dynamic simulations obta- 
ined two steady state solutions. One of them stable and the 

other unstable, by varying certain parameters in their model. 
1 5 

Kubicek et al. found in a certain range of parameters, 
three steady state solutions, two of which are stable. 

Multiplicity of solutions arises due to the depen- 
dency of certain parameters like and on composition 
and concentration of the substrate and temperature. The 
sensitivity of the stationary solutions to the individual 
values of parameters in certain cases is high. 

The dependences of the stationary values of subs- 
trate concentration on volatile acid dissociation constant 
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Table 3*2 


Values 

of Parameters 

Used in Andrews and Graef Model 

Parameter 

Value 


Unit 

Meaning 

F 

1.0 


1/day 

Flow rate 

K 

0.024593 


mol/l-atm 

Henry's law constant 


0.316 X 

10“"^ 

mo 1/1 

Acid dissociation constant 

K. 

1 

0.567 X 

10“^ 

mol/l 

Inhibition constant 


100 


-I 

day ' 

overall gas transfer 
film coefficient 

^1 

0.1 X 10 

-5 

mol/l 

Ionisation constant of 
the bicarbonate 


0.303 X 

10'"^ 

mol/l 

Saturation constant 

Sy 

0.02 


- 

Organism yield 


47.0 


— 

GO^ yield 

Sm 

V 

47.0 


— 

GH^ yield 

10;0 


1 

Digester volume 


2.0 


1 

Volume of gas space in 
the digester 


0.4 


day*"^ 

Maximum specific growth 
rate 

P 

0.93421 


atm 

Total pressure 


2.0 


^ “1 
day 

Toxicity influence rate 

Sp 

25 


1/raol 

Gonversion factor 

Gy 

Ap 

0. 001 


mol/l 

Inlet concentration of 
organism 

G© 

Sq 

04 09 


mol/l 

Inlet concentration of 
substrate 

o 

a 

o 

0.009 


mol/l 

Inlet concentration of 
dissolved GO^ 

% 

. 0;1 


mol/l 

Inlet concentration of 
net cation 

■ C»p', 

0.0 


mol/l 

Tn 1 o + r n n r on + + 1 nn o-F 
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K and the stationary valves of substrate and biomass concen- 

a 

tration on inhibition constant are shown in Appendix A. 

The stationary or steady state solutions for the 
anaerobic digester were obtained numerically from the model 
representing the digester, by setting the derivatives in the 
model to zero. 

The dependence of the stationary value of the subs- 
trate concentration on the value of flow rate, F for a value 

“5 

of = 1 0 is shown in Figure 3. 2i Three stationary solutions 
v;ore obtained. Of these three steady states, two of them 
(a and c) are found to be stable and the middle steady state 
(b) is found to be unstable. The reason is if we give a 
perturbation in any one of the manipulated variables (inlet 
organism concentration Gy , inlet substrate concohtration C* 
and inlet dissolved carbon dioxide concentration ), it 
results in the shifting of the steady state ‘b' to the stable 
states (a and c), according to the direction of the pertur- 
bations. 

Again, of the two stable steady states, the steady 
state corresponding to low substrate conversion (a) is found 
to be highly productive whereas steady state (c) corresponding 
to a high substrate conversion, is found to be unproductive. 
The reason for this is shown in Figure 3.gi The figure indi- 
cates that the rate of methane formation is very high in the 
steady state ‘cl* ^^han compared to the steady states ’b' and 
'c' and hence the best operable steady state. 




F (liters /day) 


Fig. 3.2. Dependence of Stationary Solutions on the 
of F. Klq =30, Kj=10“^, Other Parameters 
as in Table 3.2. 





Dependence of Rate of Methane Production, 
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3.4 Digester Failure 

Many cases of digester failure arise because of the 
operation not recognising the signs of impending difficulty. 

If the process variables are monitored well, it is likely 
that the failure can be thwarted. 

In a digester, failure is due to the instability 

which results from any material or effect, that interferes 

with methane formation. Three potential causes of digester 

instability and ultimate failure are hydraulic, organic and 

1 6 

toxic overloading. 

Hydraulic overloading can occur whenever there is a 
washout of the microorganisms. This increase in flov/ rate ’F' 
happens when sludge production exceeds design capacity or 
whenever too much sludge is pumped into the digester. The 
residence time can also be decreased if top deposits of 
grease and scum accumulate within the digester or when insuff- 
icient mixing produces stagnant zones in the reactor. 

Organic overloading can occur whenever there is an 
accumulation of volatile acids in the digester, that can 
inhibit the methanogenic organisms. This can result whenever 
there is a sudden increase in organic solid feed concentra- 
tion caused due to excessive digester loading (increase in 
flow rate *F') or attempting to start up a digester too 
rapidly. Another source for failure due to organic over- 
loading is when a carbohydrate predominant sludge is added. 

The decomposition of carbohydrates, in contrast to normal 
sludge degradation, in which ammonia is a metabolic product, 



29 


yields short chain fatty acids without the accompanying 
buffer action provided by ammonia. This results in the fatty 
acids remaining in the unionised form v;hich are inhibitory. 

Toxic overloading can occur when a material that can 
kill the methanogenic microorganisms is introduced. The 
materials that have caused failure includes heavy metals, 
detergents organic chemicals, ammonia and various cations. 
Fortunately, this failure can be prevented as there is always 
sufficient dilution capacity to dilute these toxic agents 
occasionally found in domestic and industrial wastewaters. 

Besides these three potential causes of digester 
instability, failure can also occur due to the effect of 
timelag in the specific growth function. The microorganisms 
present in the sludge do not react immediately to the change 
in conditions of the environment i.e. concentration of the 
substrate. The increase of inhibition by the substrate can 
also be described as an increase of the timelag (delay) in 
growth rate function. The timelag has a destabilising effect. 
For a sufficiently large timelag (tj^> 2.5 days) this effect 
is pronounced. Fortunately, in practise, the observed timelag 
is very small^Ct^^ = 0.3 day) and does not pose any serious 
problem of destabilising the digester (Appendix A). 
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Chapter 4 

RELATIVE GAIN ANALYSIS 

The traditional industrial control strategy for 
multivariable control problems, v/here several process vari- 
ables are to be controlled and several variables can be mani- 
pulated, is to use a multiloop control scheme consisting of 
several conventional controllers. In designing a multiloop 
control system, a key design decision is to determine the 
proper pairing of controlled and manipulated variables. If 
an incorrect pairing is used, the resulting control system 
might perform very pooly or even be inoperable the depending 
on the level of interaction between the control loops. 

To illustrate the multiloop control strategy and 
interaction of control loops, consider a process, which has 
two controlled variables and C^, and two manipulated 
variables m^ and m^. There are two possible control loop 
configurations. One approach would be to control C<j by adjus- 
ting and to control C 2 by adjusting m^* This is I'eferred 
to as a 1 -1 /2-2 configuration. The alternative approach is 
to control by adjusting and to control by adjusting 
i.e. a 1-2/2-1 configuration. 

Let us take 1 -1/2-2 configuration. To understand 
the nature of interaction between two control loops, the 
effect of input changes or the outputs is studied when 
(i) One loop is closed and the other is open, and 
(ii) Both loops are closed. 
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In the former case (Figure is clear that any 

change in the set point R, sp will affect not only the beha- 
viour of controlled output but also the uncontrolled out- 
put 

In the later case (Figure4ilb}, consider a change in 
set point , sp only and keep the set point R 2 , sp the same. 
Then controller of loop 1 will change ra^ to bring to the 
set point value. This is direct effect of m^ on through 
loop 1 . 

The control action of will not only bring to 
set point but also disturb from steady state. To make 
constant, the controller of loop 2 changes appropriately the 
value of m^* But a change in in turn affects output C^. 
This is an indirect effect of m^ on . 

This constitutes the essence of effect of interaction 
between two control loops. For better controls, interaction 
between loops should be minimum. 

The most prominent approach for characterising 

process interactions is the Relative Gain Array (RGA) method, 

23 

proposed by Bristol. This method, based on the steady state 
process model, provides useful information to categorize 
process interaction. However, process dynamics needs to be 
considered in many situations and another method, dynamic 

OR 00 

gain array method was proposed, * which gives a useful 
insight into the dynamic behaviour and also provides a 
measure of dynamic interactions as well. By studying both 
the steady state and dynamic approaches, a total picture of 
the behaviour of the process can be got and measured. 
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4. 1 Steady State 

To illustrate the relative gain array method, 
consider the steady state process model 


'SS 


= K 


^S 


(21 ) 


where Cgg and m^g are the steady state changes in the n 

controlled and n manipulated variables, respectively, and j< 

is the n X n steady state process gain matrix. 

"til 

The relative gain between the i^ controlled vari- 
"bli 

able, and j manipulated variable, m. will be denoted 

3 

This dimensionless quantity is defined as 


by Ny. 


. ^en loop gain 

^ij Closed lo'op gain 


Mathematically, it can be expressed as 


( 22 ) 




(23) 


In this equation, the partial derivative in the 
numerator is evaluated by keeping all the manipulated vari- 
ables except m. constant (Figure 4« la)* Similarly, to evaluate 
the partial derivative in the denominator, all of the cont- 
rolled variables except are held constant (Figure 4.1b). 
This later condition can be achieved (in principle) by adjus- 
ting the other n-1 manipulated variables. 

From eq. (l ) it is apparent that the open loop gain 

between C. and m. is K. an element of the process gain 
X . 3 _ 


ma 


trix, K. Bristo^^as shown that the closed loop gain is 
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1 j where is an element of the matrix inverse of the 
transpose of K. 


<s» r » *1 

K = K. . = [K^3 


Now, the expression for N. becomes 

xJ 


Tn-I 


(24) 




K. . K. , 

3-3 3 3- 


(25) 


The Relative Gain Array is defined to be the n x n matrix. 


A 


[^..3 

3.J 


The relative gain array is a square matrix which implies that 

number of manipulated variables is equal to number of contro- 
04 

lied variables. 

A useful characteristic of the relative gain matrix 
is the fact that the sum of the of any row or column of 

the matrix has a sum equal to the unity. 


N 

•^1 ^3 
1=1 


1 


for j = 1 , 2, . . . ,N 
(summation by rows) 




= 1 


for i = 1 , 2, . . . ,N 

(summation by columns ) 


N 
2 

•A 

3=1 ■ 

The relative gain gives an useful measure of interaction and 
provides a recommendation concerning the proper control confi- 
guration, In particular^ 

(l) If Xj - = 0, then manipulation of m. has no effect on the 
control variable C^. 
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(2) If 0 < . <1, then an interaction exists and manipu- 

lation of variables other than m. also affect the steady 
state value of the control variable The smaller the 

value of the larger the interaction. 

(3) If = 1 , then the control loop formed between C. and 

1 

m^ does not interact vdth the other loops formed between 
the other controlled and manipulated variables. So 
manipulating other variables (other than m ^ ) does not 
affect C^. _ 

(4) If . < 0, then other manipulations (other than m.) 
causes a strong effect on and in a direction opposite 
from that caused by manipulating m^. The interaction 
effect is very dangerous in this case. 

(5) If A.jj > 1, then the response of the output C. is held 
back by the interaction from the other loops formed 
between the other controlled and manipulated variables. 
The longer the values of the relative gain above unity, 
the larger will be the holding back effect. 

Criteria for Selection of Control Loops* 

From the characteristics of relative gains and rela- 
tive gain array, it is seen that the best control loops (with 
minimum interaction) are got by pairing the controlled outputs 
with the manipulated variables m^ in such a way that the 
relative gains are positive and as close as possible to 
unity. 



4*2 Dynamic Approach 

Although the steady state RGA method provides useful 

information, it neglects the dynamics of the process. In some 

situations, the dynamic characteristics may be an important 

factor in determining the control configurations. Consequently 

other measures have been proposed which includes both static 

25-^27 

and d^mamic behaviour of the process. 

We have adopted the Witcher and McAvoy approach. 

Witcher and McAvoy’ s measure of process interaction was based 
on open loop step responses. This information can be got 
from experimental data or from transient response of a process 
model, for e.g., a transfer function or state space model. 

Thus this method is widely used as it does not require a part- 
icular type of process model. 

The interaction measure is based on a ’dynamic poten- 
tial’ which is the integral of the open loop step response 

C. (t) to a step change in m. at time t = 0, 

4- J 

0 

0..(e) = / C.(t)dt (i,a = 1, 2,...,n) (26) 

^ Q ^ 

Here the selection of 0, the time period over which the 
integration is performed, is somewhat arbitrary. It is speci- 
fied as 20^ to 100^ of the dominant time constant of the 
process. By substituting '0. . (©) for steady is tate-gain K, . in 
equation (25), Witcher and McAvoy constructed a ’Relative 
Dynamic Array’ (RDA) with element X. .(©) given by 
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where is a element of 0(9) and 0(9) = 0(9) *^ 

The selection of loops in this case is similar to 
RGA method. However, this method has a shortcoming in that 
both \^j(9) and the recommended pairing depend on the value 
of 9, which is not well defined. 

A new measure of process interactions, a modified 
form of Witcher and McAvoy method, was proposed by Gagnepain 
and S eborg'. 

This v;as an extension of the RDA concept since it was 
also based on the open loop step responses. However, to avoid 
the disadvantages of the RDA method, the resulting dynamic 
response data are used in a different manner. 

Suppose that the process is initially at steady state 
and that a unit step change m^ occurs at t = 0. During the 
time interval (O, d. j), C. is not affected by m. and the 
'dynamic gain' of the process at this time period is zero. 

For the time interval (d^j, ©), the average dynamic gains, 

D*., between G. and m. can be calculated as 

„ (Average change in C. ) 

ij “ * (Change in m. ) 

, 1 


But since, at steady state a unit step change was made in the 
process, it follows that Dt. is equal to the average of C. 

Ij 1 

over the time interval e). 




9 

/ C.(t)dt 

j 1 


k pi ecewis e cons tant f unction is 


defined conveniently as 
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D. ^(t) = 0 


= D*. 


for t < d. . 

\J 

for d. . < t ^ 0 

J- J 


In a further modification of Witcher and McAvoy method for a 
process with negligible timelag, we have the average dynamic 
gain, D. .as 


D... / G. (t)dt 

^3 0 - ^ 

a 

where © is the largest time constant in the process transfer 
function matrix. 

In analogy with RGA and RDA, average dynamic relative 
gain is defined in terms of D. .(t) 

p..(t) = D..(t)D..(t) 

where D.^ is an element of D = 

sj 




39 


Chapter 5 

RESULTS AND DISCUSSION 

5. 1 Ca>libration of the Simulation Software 

In order to check the consistency of the software 
developed we tried to reproduce the batch simulation results 
of Andrews and Graef using his data. Tremendous difficulty 
was encountered in producing this match. During this diffi- 
cult period several parallel numerical techniques were tried* 
Powell algorithm, Brown's method and Math proton method. 
Continuous debugging of the program led to discovery of an 
error in Kubicek's model presentation. Dynamic equation 
for partial pressure of carbon dioxide carried a plus sign 
with the pressure term, whereas the physics of the problem 
indicates a minus sign. 

After removal of the error indicated above, dynamic 
simulation was tested and confirmed with ^ubicek’s simulated 
results (Appendix B). 

5.2 Steady State Analysis 
5.2.1 Hydraulic Loading 

From the equations representing the model we have 

Cv » C* » as the inlet concentrations (manipulated vari- 
■^0 *0 ^0 

ables), and Cy^, Gg, C^, P^, C^ and C^ as outlet variables. 
Among the outlet variables, the controlled variables are to 
be chosen. In the model, equations (19) and (20) have a 
linear first order form and C^Ct) and C^{t) are stable and 
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and tend to the inlet conditions Cy and G™ • Thus, accor- 
ding to the model, only the first four variables Cg, 
and pQ need to be controlled. 

Of these four variables, Cg, and are measure- 
able and also can either be measured or microorganism 
activity can be estimated from the methane gas production 
rate. Control strategies should be developed using relative 
analysis such that each manipulated variable controls one 
controlled variable only so as to achieve a near perfect 
decoupling of the loops (with minimum interaction) formed bet- 
ween the controlled and the manipulated variable. 

The current practise of digester operation involves 
empirical control rules using rate of methane production QCH 

/j. 

and pH as feedback control variables. Considering them among 
the variables to be controlled, we. have 6 possible process 
variables to be controlled, C^^, Cg, C^ and P^^ from the equa- 
tions representing the model and QCH^ and pH from current 
operational practise. 

As mentioned earlier, flow rate 'F’ is considered as 
a disturbance variable. Similarly Cy (which causes toxic 
overloading) is also taken as a disturbance variable. 

So, a control strategy should be made which holds 
good for various values of the disturbances. 

Such a control strategy can be realised by perf arming 
th« relative gain > analysis for the system. 

Similarly, relative gains were computed for a parti- 
cular flow rate (Table 5.1) by making a 10?^ perturbation in each 
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T A 

b L E 5.1 


RELATIVE 

GAIN 

array FJR 

the 20, CniJTRnL 

SF.TS ObTATSiRO 

AT STEADY STAIE 'd' FOR 

A FLO^‘ RATE Ot 

F=1.U LTiERS/OMY 



cxo 

• 1 ■ 

cso 

CCO 

(1) 

CC 

PC 

PH 

-14.90452 

13,R6fi03 

2.636417 

-109,0395 

1 1 1 .0603 
-1 .02999 

124.9436 

-123,9374 

-0.06426 



CXO 

c:su 

CCO 

(2) 

CS 

CC 

PC 

0,70903 
-1 1 ,6048 
12.09580 

0, '2 89870 
-112.1218 
112.63200 

0,001090 

124.9276 

-123.9278 



CXO 

CSO 

CCO 

(3) 

CX 

CC 

PC 

1.115840 
18,47280 
-1 8.5RB6 

-•0.160200 

-142.2297 

143,34570 

0,0001 80 0 
124,75680 
-123.7570 



CXO 

CSO 

CCO 

(4) 

CC 

PC 

OCH4 

11231.064 

13.868030 

71.430209 

-11291,992 

1 11 .069300 
-'70.539901 

124.94380 

-123.9374 

0.1097102 



CXO 

CSO 

CCO 

(5) 

CC 

0CH4 

pH 

-1,12146 

0.087545 

2.033922 

1.343880 

0.69627? 

-1.040158 

0.7775b 

^.216178 

U.U06236 



CXO 

CSO 

CCO 

(6) 

CX 

CC 

PH 

0,47677 

-5.543077 

1.166307 

0.399252 

5.1 75179 

-• 4.57433 

0,123977 

-^.532101 

4.4061200 



CXO 

CSO 

CCO 

(7) 

CS 

CC 

pH 

5.548716 

9,351408 

-13.90012 

-18.27044 
85.218750 
-6 5. 94 830 

13,72172 

-95.57015 

8u, 848430 
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cxo 

cso 

ccc 

(a) 

PC 

OCH 4 

PH 

-1.12839 

0,094661 

2.033710 

1.352230 

0.637301 

-1.040534 

O.776i501 

0 . 2 1 7 b 3 6 
0.00631 54 



CXO 

C50 

CCO 

(9) 

cx 

PC 

pa 

0,49y269 
—0 ,625 34 
1.127070 

0.33U51 

5.032660 

•4,41383 

0.1 20569 
~i,40 / 3-2 
4,2b675o 



CXO 

CSO 

cco 

( 10 ) 

cs 

PC 

pH 

3.409475 

5.346540 

-7.75604 

-10.25540 

48.721860 

-'37,46940 

7.81593 i 
-5 3.071 ‘.3 

4 6,2254 00 



CXO 

CSO 

CCO 

( 11 ) 

Cs 

0CH4 

Pri 

0.594170 

a,u781b8 

0.327658 

0.2927301 
0.7074331 
- '0 , '0 0 0 1 / I 

0. 1 13u9 

0.?14:>9 / 2 
0.67251 24 



CXO 

CSO 

CCO 

( 12 ) 

CX 

oca 4 

pH 

1,117674 

-5,11758 

0,000006 

-'0,140041 U 

0. 94D5D601 

0 . 1 9 9 5 3 0 0 1 

0,0223601 
C. 177 1760 
0.8004566 



CxO 

CSO 

CCO 

(13) 

cx 

cs 

pti 

0.44dQ93 

0,357022 

0,l968t»0 

0.4 2506001 
l.lsl256Cl 
- 0 . D 0 S 3 3 0 2 

0 , 1 2 8 b 4 U 1 
-U-53d23^ 
1 . U' » -i 4 j \ 



CXO 

CSO 

CCO 

(14) 

cx 

PC 

0CH4 

1.117678 
0.(/00003 
-0. 1 17d8 

-0,1 175000 

0. 51767201 

D.'8 9932B0l 

«u. OOul n 
0.^623241 
0.2176530 



CXO 

CSO 

CCO 

(15) 

CX 

cs 

PC 

0.435050 

0.432591 

0.132360 

0 . 432 o 500 i 

1 .36938000 
-'0. 3014 37 

0.1i2S9oC- 

-OeOiO/O 

1 .^d907 70 
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CaO 

cso 

r» --o p 

C16) 

CX 

CS 

CC 

0. 439863 
0.429521 
0.130605 

0.4290001 

1 .351 7660 
-■0 . 7 90 77 5 

0. 131177 
-0.79179 

1 . foOl?! 



CXO 

CSO 

n - 

(17) 

CX 

CS 

QCH4 

1.317690 

-0.00001 

-O.3l7o8 

-O.OOJll? 

0.4924R01 

0.5076330 

-0. 317570 
0.5C7525 
0. R10051 



CXO 

250 

cco 

(IB) 

CX 

CC 

0CH4 

1.11767B 

O.O0u004 

-0.11768 

-0.117502 

0.2153050 

0.9011960 

-U.ODul 75 
0.7 o3o971 

0.2164830 



CXO 

CSO 

f'* n 

w- w- V' 

(19) 

CS 

PC 

0CH4 

0.60B273 

0.215358 

0.176368 

0.3927850 
-0, 000221 
0.6074351 

0 - u 0 i i ) 3 8 

0 , 7 B 4 0 ^ 0 


CXD 

CS 0,b3b286 
CC 0.^16712 
QCH4 0,175000 


cso 

0 . 2 9 ? 7 8 b 1 
»0.oa02/:2 
0.^074870 


rrf'i 


-0.00x077 

0,7685091 

0.2175025 


C203 
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of the three manipulated variables Cy , and . This 

■^0 *0 ^0 

yielded 60^ 20 possible relative gains. The best confi- 

guration should be chosen from these twenty possible relative 
gains. 

The following criteria for discarding certain control 
loop possibilities was used. It can be seen from Table 5.1 
certain controlled variables, when they exist together in a 
relative gain array or a control configuration yielded control 
loops with high interaction between the loops, as denoted by 
their relative gain values. They should be discarded. 

The following sets of controlled variables, if exis- 
ting together, were found to give high interacting loops 
(see Table 5;1 )i 

(a) Cq, Pq and pH as seen in Set No. 1. 

(b) Cq and Pq as two of the three variables in a set. This 
occurs in three cases (Sot No, 2, 3 and 4) 

(c) Cq and pH as two of the three variables in a set. They 
occur in three cases (Set No. 5, 6 and 7) 

(d) Pq and pH as two of the three variables in a set. They 
occur in three cases (Set No, 8, 9 and 10). 

Thus from the twenty possible control configuration, 

10 are discarded this way. The reason for this high inter- 
action is that Cq is the concentration of dissolved carbon 
dioxide, Pq is the partial pressure of carbon dioxide and pH, 
depends on the concentration of carbon dioxide. This close 
relationship between them make them yield control loops with 
high interaction. 
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From the remaining 10 relative gain configurations, 
the best ones are to be chosen. For this, the best variable 
among P^, pH and should be considered, as one of the vari- 
ables of the relative gain array. These three variables are 
interrelated and behave in a similar manner. Therefore it is 
difficult to choose among the three. 

pH is considered ds the best measurable among them, 

Pq is also easily measurable and is comparatively difficult 
to measure. This reduces to three controlled variable sets 
as seen from TableCS-O Chayirvi^ 

(i) QCH^ and pH, Set No. 12 

(ii) Cg, QCH^ and pH, Set No. 11 

(iii) Cj^, Cg and pH, Set No. 13 , 

In addition to these three sets, two other sets having 
P^, i.e., (C^, QCH^, P^) and (Cj^, Cg, Pg) are also considered. 

The other remaining five sets have been rejected. The 
rejection lies either in the set having a variable which is 
comparatively difficult to measure or in the form of a poor 
relative gain array, denoting high interactions between the 
loops. 

The relative gain analysis is repeated for different 
flow rates, each time yielding twenty relative gains^. All 
the configurations got with each of the flow rates were similar 
in nature. Using the same criteria for discarding control 
loops, the five control configuration selection was confirmed, 
for all flow rates. 
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In order to select the best among these five control 
configurations, an index is used formulated as 


L.S.N.I. 


I (X? - 1)^ 

i=1 


We called it the least-square non-interaction index where 'h. 

e- 

is the best possible relative gain got considering that only 
one manipulated variable controls one controlled variable; 

The smaller the value of least square non-interaction 
index, the better is the control configuration. 

Figure 5.1 is a plot of least square non-interaction 
index versus flow rate, F. The plot indicates that there is 
a sharp change in least square non-interaction index for each 
of the configurations around a value of F = 1.35 lit/day. 

We have already seen that the rate of methane production 
increases and reaches a maximum value at F = 1.5 lit/day 
(Figure 3.2). To exploit this methane production to the fullest 
extent, pairs of control configurations are considered where 
in a configuration is used for a range of F from 0 to 1.35 
lit/day and its pair is used for a range of F beyond 1.35 
lit/day. The 'configuration pair’ which operates with least 
square non-interaction index would obviously be the best pair 
of configuration available. 

This is easily seen from Figures 5.2and5;3. Three 
sets are considered. 

(l) Split range control i Configuration /QCH^-Cg / 

pH - used for a range of F from 0-1,35 lit/day, 

^0 



10.0 


0.01 


0.00 lL 
0.2 



0.1 — P 

/ 

- / 

_/ 


0.6 0.8 1.0 1.2 1.4 

Flow Rate, F (liters /day) 


Fig. 5.1. Plot of Least Square Non Interaction Index 
of Gontrol Configurations for Different 
Rates . 


Least Square Non Interaction Index 



Fig. 5.2. Plot of Least Square Non Interaction Index vs 
Flow Rate, F. 
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Spilt Control : 
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Cx~Cxo^s~Cso /pH-Cco 
(beyon(in.35 l/d) 

Full Range Control ; 
Cx"Cxo ^Q.ch4"Cso/Pc“Ccc 


0.001 


0.4 0.8 1.2 1.6 2.0 

Flow Rate, F (liters /day) 

Plot of Least Square Non Irfieraction Index vs 
Flow Rate, F . 
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sequenced by configuration - Cv /Ce “ C*. /P„ - C- 

A ^0 ^ *0 ^ ^0 

operating for a range of F beyond 1.35 lit/da y, 

(2) Split range control? Configuration Cv - /QCH. - C^ / 

X Xq 4 Sq 

pH - Cp operating for a range of F from 0 to 1.35 lit/ 
day sequenced by configuration Cv “ Cv /C* - C^ /pH - Cp 

X -^0 ^ *0 ^0 

operating for a range of F beyond 1.35 lit/day. 

(3) Full range control? We have a single configuration 

-* Cj^^/QCH^ - Cg - Cq operating in the full 


range of F (O to 2 lit/day). 

The Figures 5.2 and 5.3 indicate that all the three 
sets can be employed with equal benefit in controlling the 
anaerobic digestion process. 

However, controlled variable set Cj^ - ^^^4" seems 
to provide the best control since it is operable in the whole 
range of F and it also eliminates the difficulty of changing two 
controlled variables (other than Cj^) beyond a particular 
range of the flow rate. 


5 . 2. 2 Toxic Loading 

The relative gain analysis was performed for different 
values of the second disturbance variable, C^ • As before, 
the five control configurations were confirmed and compared 
using the least square non-interaction index (Table 5,2)* 

Figure 5i4 indicates that the toxicity has virtually 
no effect on the configurations. 







Table 5.2 

Least Square Non-Interaction Index of the Five Configurations 
for Different Toxicities and Flow Rates 


Flow rate F 
(l/d) 

Toxicity, C~ x 
0 

10“^ (m/s) 



0 

1 

2 

3 


Cr/ — Cy /C^J ~ 

-•V Aq O 

JSq ^ 

^0 


0.6 

0.8 

1.0 

1.2 

2. 39183 
1.52991 
0.90327 
0.48337 

2. 36981 
1.49780 
0.88497 
0.47319 

2. 34624 
1.47097 
0.86706 
0.46325 

2. 3221 0 

1 . 44274 
0.84933 
0.45362 


Cg - C^^/QCH^ 

- C« /pH 
*0 

■ 

^0 


0.6 

0.8 

1.0 

1.2 

0.19389 

0.27320 

0,35751 

0.45691 

0. 1 9644 
0.27590 
0.36055 
0,46097 

0.19912 ■ 

0. 27864 

0. 36363 
0.46513 

0.20103 
0. 29843 
0. 36677 
0. 46943 




- Cp 
-"0 


0.6 

0.8 

1.0 

1.2 

0.04015 

0.05002 

0. 05720 

0. 06677 

0. 04067 
0.05057 

0. 05774 
0.05821 

0. 0441 1 

0. 05102 

0. 05836 

0. 06982 

0.04160 
0, 051 96 
0.05903 
0.07161 


Cy *“ /Cf ~ 

A ..Q ... 

Ce /pH - C 
*0 



0.6 

0..8' 

1.0 

1.2 

1 . 44082 

0. 90524 
0.50731 

0. 30279 

1.42961 

0.88850 

0.49965 

0.30390 

’ 1 . 40638 
0.87217 
0.49240 
0,30412 

1.38262 
0.74528 
0.48389 
0. 30426 



“ Ojg /Pp ■“ Cp 
^0 ^ ^0 


0;6 

0.8 

1.0 

1.2 

0.04550 

0. 05390 
0.06123 
0.09128 

0.04613 

0. 05989 
0.07187 
0.08147 

0.04684 

0.06049 

0.07255 

0.08671 

0.04771 
0. 06607 
0.07326 
0,08812 


Least Square Non Interaction Index 


1.0 
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Fig. 5.4. Plot of Least Square Non Interaction Inde 
vs Toxicity, Cjq at a Flow Rate , F =1.0 
liters /day . 
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5.3 Dynamic Analysis 

The dynamic relative gain analysis was performed 
using the dynamic model of the digester. This was done by 
simulating the response of the controlled variable to the 
change in manipulations and then computing the dynamic rela- 
tive gain by using the modified form of Witcher and McAvoy 
method, described earlier. 

Figures 5.5 to 5. 9 represents the dynamic responses of 
the six controlled variables to the change in the manipulated 
variables. 

Table 5.3 compares the least square non-interaction 
index of the five control configurations for both the steady 
state and dynamic analysis. It can be seen that the least 
square non-interaction dynamic index is a shade better than 
the least square non-interaction sta.tic index. Also, it 

confiims that the ~ ” ^^0^ 

QCH. - Cc /P^ - C„ configurations are the best possible for 

an anaerobic digester. 
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Table 5;3 


Comparison^of Least Square Non-Interaction Index for the Control 
Configurations in Steady State and Dynamic Analysis 


Control configui'ation 


Least square 
non-interaction 
index in steady 
state analysis 

oiV r-x. VC? V-Vs/Vs*v 


Least square 
non-interaction 
index in 

dynamic analysis 

o-t P \ • o It 



Cy /QCH. - Ce /pH 
Xq 4 Sq 



0.35751 


0. 26048 



a^/QCH^ 


Ce /pH 
*0 



0.05720 


0. 02925 


Cy - Cy /C« - Ce /PH - Cp 
X ^0 ® '^0 


0. 50731 


0. 331 27 


” '^Cq 0.90327 

/QCH^ - Cg /P^ - 0 


0.56231 


0. 061 23 


0.03085 
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CONCLUSION 

System engineering tools like the steady state and 
dynamic relative gains succeed in their attempt in quanti- 
fying the process operation and developing a more rational 
control strategy for the anaerobic digestion process. 

The results of the relative gain analysis bring out 
the following facts: 

(1 ) The present day control practice is confirmed: 

(a) The present day operation is in controlling the pH 

by scrubbing of carbon dioxide and recycling it back to the 

digesteriE-elative gain analysis gives a configuration 

Cv - Cv /QCH. - Cq /pH - Cr. in which pH is controlled by 
X Xq 4 i>Q Cq 

manipulating the inlet carbon dioxide concentration. 

(b) Another recent practice is in controlling the rate 
of methane production by recycling the concentrated sludge. 
The relative gain analysis suggests that the rate of methane 
production is to be controlled by manipulating the inlet 
substrate concentration, which can be done by recycling the 
concentrated sludge. 

( 2 ) New controlled possibilities are suggested: 

controlling pH by 

scrubbing and recycling of carbon dioxide does not hold 
good in situations where hydraulic loadings exceed a parti- 
cular value. A new controlled possibility - Cj^^/QCH^ - 

— wherein the partial pressure of carbon dioxide 
^0 
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is controlled by manipulating the inlet carbon dioxide concen- 
trations holds good both for hydraulic and toxic loadings. 

The same benefit can also be achieved by having split control: 

(a) Configuration "* ^Cq 

ting upto a hydraulic load of 1«35 1 it/day sequenced by con- 
figuration % - /C 3 - Cg /pH - G^^operating at hydraulic 

loads beyond t*^ lit/day. -It should be noted that is a 
change in the controlled variable from QCH^ to Cg during the 
sequential operation. 

(b) Configuration - Cj^^/QCH^ - Cg^/pH - opera- 

ting upto a hydraulic load of 1 ,35 lit/day sequenced by con- 
figuration Cjj - - %q/^g “ ^Cq beyond 

1.35 lit/day. Here, it should be noted that two controlled 
variables change (QCH^ to Cg and pH to Pq) during the sequ- 


ential change. 
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SUGGESTIONS FOR FUTURE WORK 

The following are the suggestions for future works 

\a> 

and Sv should be identified. There exists an uncei>- 
tainty in the substrate composition and the associated 
parameters and in order to have faith in the control 
parameters, these process parameters need to be identi- 
fied by employing the input-output data of the process. 

The relative gain analysis of the anaerobic digestion 
process aids the multivariable control pairing decision, 
however, other control strategies like the adaptive control 
and feed forward control can also be tried for the anae- 
robic digestion process. 

The relative gain analysis has been done for the conti- 
nuous feeding of the digester (represented as a CSTR), 

The analysis also needs to be done for batch feeding 
and periodic feeding of the digester. 

The current model of the digester considers acetic acid 

22 

as the lone substrate. Recent research indicates that 
other higher fatty acids can a'lso produce more methane 
gas by providing the proper biological environment. So 
focus should be done. on modelling and simulation using 
different substrates, individually. Not much data is 
available in the literature and this could be a great 
helo for the operation of the digester. As a corollary, 



The process parameters like 
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modelling and simulation should also be done for diff-*- 
erent substrates coexisting together* 

(5) Membrane s^eparation process like dialysis can separate 
the two phases — the acetogenic phase and the methano- 
genic phase, co-existing in a digester. If model is 
available for both of these phases, a combined control 
strategy for both of these phases can be developed. The 
model for the methanogenic phase is currently in use 
since it is considered as the rate limiting step of the 
digestion process. Our focus should be on modelling the 
acetogenic phase. 
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APPriMDU d.l 


PRDGPA'I F 3 R GPTj'TMG SrM’lJiJAPV snb'U'IjNb 
By SOLyiMG NOM-tlNFAR A’T’I 


THF PARaMRTcRs USciD are 


FC = FI.DW RATE.L/DAy. 

AK = riFr'»R/ S bA« C^NS^^.^^T, MJL/r. AT'*. 

Aka=actd uissbri Arr j'j :jN!brA'rr,r.‘n,,/i,. 

A^I=:I,NHi3ITl^,^. CnuSl’ANT, MJL/r,. 

AKLA = Dyfc:RlvI.b GAS rRAMSPbR FILM Cni'.F j' , /n aY . 

AKl =iaMiSATinNi c^NS^A^l.^ of bicAKPjWATt, - iDg/u.- 
A^S = SATURAT10N S f A \ T , '.v)t,/I. . 

SXS=ORGAMlSvv yiFLD. 

SCX=CD? yiELP 
SMX=Ch4 yield. 

VsDIGcISTER bljUjD VJfiJME^L. 

VG^VDLID'IE of gas SPA:E,b, 

HMU=MAX. specific GHHaTH RAi'P-r^/uAy. 

ARTsTOXIClTY INFf.UEvCE HAi’E,/uAy. 

SP=CnNVERSl3F FACr3R» b/MTL. 

P=TaTAl. PRESSURE, atm. 

CXO,CSO,:CO,C^O,CT^=I^iLET CUMCF.i,TKAj‘Ii)+;S IIF DnG J. S'lS , M ‘U TJ 

substkate, dissolved co2,?jet CATin.; A'’u iDXTC Rc-sr-ctj. 

niMENSlO'J Xl30) 

COMMOF FOtAK,AKA, AKT, AKLA,AK1 , aKs , SXS , SC X , S'U , w , vG.Hv.-J, AXT, St-', 
PCXOCSDc''OCZO,'’rD 

6 aTA'fU,Ak 7 A^A,AKI, AXr.A, AKl ,AKS,SXS,SCx,S’'U, W , v r, , h >. 11 , ft K T , S R , p , 
CXO,CSO,CCO,CZO,CrO/.D.6,0.024b9 3,u. J1 bF-4-, 1 . U K-j , u 0 - D , 3 . 1 E-d, 
0,33 3E“4,0,02,47.0,4 7,D,10.0,?.€, 0.<; ,?.0,75.t>,u.^3^7l , , 

0.09,0.0081 ,0.1,1. 

0PFN(UNir=3l , DEVI CE=C USE* 1 
X(l)=:0.002b 
X(2)=0.010 
X(i) =0.001 
X(4)=0.0i9 
TYPE^, (X(I) ,1=1 ,5) 

rj=:4 ; F UMS T G = fa ; M AX T r =2 5 ; T PR I M r = l ; E ps= 1 . E-R ; V D = 3 1 

CALL NUMLlwCf*, mOMSIG, viAXlT,IPFl'.r,X,Ers, ASibn, .DJ 


STOP 
END 

subroutine 


M DN Ll N ( N , N IIMSIG , MAX I T , 7 pR 1 F l , X , EPS , aS i R k , s, !) j 


This siioroutlne soWes a system ot a sinuitanioJK n^n l 1 near p 4UC, 
tions.Tae metriod used i s at least auadrat ica-i i y con« er^enL aTJ 
requires only (n**2 /2* 3 + N/2) function e^/aiuatiois par ItefaLive 
step as Compared wltn (Mt2 + w) evaiuatioas for Ne.kun.'-'s iiecto''’. 
This results in a saviaqs of coii.putat ionai eitort far safitcient, 
complicated functicns.' The method does not r^auire t.ne user 1 0 
.furnish any derivat i-ve.s . 

INPUT PARAMETERS FDLLJ/' t 


N = Mumper eaations (= number of anKnownsi 
NUMSIG = Musiber Of significant cjlniLs desired. 
MAXIT = Maximum number of iterations to be ar.ed. 
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IPBINX = OuLPJt optlnnoutDUt Jf = 1 » ■ 

are always output ; ;*f\XlT exceeu un,i siuoular J-r.>Di:^ . 

X = Vector of prat j on wDl teruJ^acOa if 

iti tPe svsteP. 


i. - 9 ^ ^ .'4 > .-.v 


u JL ;» j. 'laj • * 4. » .A 'V- « * • - 

of toe f'Jpction values IS satistie>i. IP fur- ^ 
terminated oV one of toe criteria, simpj^ s^u 

verv stringent. 


OUTpUl PI^R^'>lF rERS FJLb^rt ; 


?^''=%ol 3 ?lon'’o£‘tS’sys?l/(or'P.‘sl- wproxi"-' «" 

^ I r" r' / '.-rn 


REAL XC 3 O),EARlCi 3 ),rF-iP( 30 ),ri)Fl. 3 b ,31 ) , RFu-''^'* • ^ 'F-^ i ' 

Ks 5 k“'j|o"Bpf^iy^?Sp 00 . 30 , 

rof^ '^nf-i/ V ar/xR » , xs , X I 


DELIA WlLu BE A FUNCm^ OF TBF m,,CHT.F A lO i P.K.-^oSl' 


TXPE»,cx(na=i» 5 ) 

RELCOM = Id.E + Ot-* t-MO^lSlC) 

Do’^TOO M il,MAXlT 
JQUIT =0 
FMAX = 0 . 

F 0 RMATn|t 7 L;lb.H/Cp 3 J.^,K^ 
LooKitpn » j 5 =o 


I.OOXUFti.u.-.. ..Ptial eftP-t «itno.t h.vin. tn 

Tne arrav LDORUP ^ partiai 

IntercPange rows or columns. 


no bOO^ *^ 7 i -iA 1 a i 
IF(K-l) I 3 ‘ifl 34 rl 3 i 

Km, riiK'kmi»,n,x,is!13, COE, LOOKUP) 


set op partial derleatlees of Kth funcllPa.. 


1345 
1 35 


CALL AJXFCM ( X. Fr KV, ■ 
FMAK= lFvl|x, AB^F ) ) 

IF(ABS(F) .LE.FPS) lU i>-i3 

I0niT = IQL!lT + l _ 

IpClOUTI.RE.f- 1 P’Ul-l 1345 

G0T0^725 nn 

FAClOR = ,001L+00 
ITALLV = 0 . 
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151 

Jt)l 


190 

200 


202 


203 


22G 

500 


r SlGivT 


209 

210 

205 


/O J 


DU 200 I=K,N 
IXFMD = liOnKOp(K,l) 

HOLD = X(1TE''.P) 

PR« TE^ME°fU»CT13S JF If 
Bt COMPOTED AS ^ 

LINS aH’W An 8 DIGIT 
f^f-rr'^FACTOR^ABSCHDLD) 

H — AMInI (F?lAX» ETA) 

IFCH.LT.PHEC) ,H = PR=-- 
XCIl’FrtP) 

PAMrBACK{KMiN^N:X,lSU3, CDF. LOOKUP) 

PARI ( 1Te5p)'^^=" CFPLJ^-F J /H 

iH APS(PARTIT IE«P) ) .or. DELTA) r iTD ' 

TFC APS(F/PART(T1EMP) ) .LE.l .L+i5) GOTO 
irALl'lf = ITALLY + 1 

?l?(i?'-t“EKLt.N-K) Ef l 2»2 

FACTOR ^ #0 E t ^ ^ - .. 

ifcfactor.gt.iu S-ITJ 77b 
GOTO 135^ P^T.-! 

H-UPS(F^^>T^Jr^:MP))!M.^FLTM go;.) ,775 

COEtK,N+l )=0.0F+00 
KmAX = IIEMP 

goto 500 

Find partial deratWBS of largpst aosolute 

Sp«AX%TBl(PARTfK««) 

r^-Pl.OE.N 

If( IEST.OT.DERMAX) -jJTJ 209 

LOOXHPCKPL^Si 1 )=Kf‘iAX 

KMAX„=^JSUB 

LOOKUP ^ ^ 

JF(ABs'(PART(KMAX)).e0.:).O) GOin 7 75 . 

^^4 ««r,»c fTT Kth row of triaadJla^ linear 

!f F values of);(l,) 

nn ?20 J = Kf^L^^S/N 

XtK.viAX) = CDFlN.N+l } , 


veil?! 


jj S" 

■E'^E DEa- 


syste^' 1 
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610 

625 

630 


649 

650 
660 
700 


1753 
C 334 

1763 

- 

C 3 35 

725 

"7777 

750 

7768 
C 336 

6777 

751 

7515 

C 337 
C 338 
775 

W 

752 
7525 

C 339 

C 401 
800 


If (N.EO.l) GQT3 610 

CALL RA:k(N“ 1 f .'J,X£l5jB,C3E,LdnKnp) 

IF(rt-l) 650,650,625 

Test for coavergenr 

DO 630 I=i,\ 

IF(ARSCTEmPCI)-X(I) j.GT.465(X(Tj) + Rh:LC0:\(') OO.p 61P 
CONTINUE 

JTEST = JTEST + 1, 

IF(JTEST-3) 650,725,725 
JTEST = 1 
DO 66u I = 1,N 
TEMp(l) = XCT) 

CONTINUE 

npF,4(UNlT=ND,DFVICE=/05K' ) 

WRlrE(NL),l753) 

FOPmM(/5x , 'NO C0NVER3E.\CE, MAXIMUM lOT'.BER OF iTr.RHTiOsiS 
ASTER=1 . 

TYPE 33^ 

• FORMM(/5x, 'N 3 CONVEPJENCE. MAXTi-.U-i RU-’.HEH OF 1 T i-.R kT I S NSFj. 
IF( IPRINT-NE.I) GOTO oDO 
WklrF.CND,l753 ) 

FDRMAKSX, 'FUMCriDN VALUES AT T,,K I.aST APPR3> t hT aO^. FOLI-j'-,: 
TYPE 335 

FDRmaHsX, 'FUNC i'TDN VALUES AT 1 r.F l.uSI’ APpR JX i p AT in FniJ.j'’ : 

IEIjAG— 1 
GOTO 7777 

IFdPRTNT.N'E.l) G3r3 900 
OPEN(UNiT=ND,DEVlCE=.'3SK' ) 

DO 750 K = 1,N 

CALL AUXFCNCX,PARrCK) 

CONIINUE 

IFdFLAG.NE-l ) GOTO 8 777 
WKTrEtf'.u,77B8) (PARnK),K = 1 , ,v ) 

FuRMATCbE20,8) 

TYPE 336, (PARTCK) , X=l , ‘J) 

FORkAT(oE20.B) 

GOTO 800 
WRT1F(N0,751 ) 

FURMAT(//'C3N’VERGEW:E has bee,. ACHTEVED.T-iE I- .iR PT i P , Au ' 
KRlrFC'NL),75l5) (PAftrCK),K = 1..) 

FDRMATCSX/aT TliE FINAL APPROaI A JT uf! F3r. uO-^ : ( :»E/. 0 . 9 J / ) 


V' AL.’’tS. ' 


astep=?. 

type 337 ^ 

F0RMAT(//'C3N VERGENZE has been achieved. THF diKCTiP 
TYPE 338, (PARK K) , !4=1 , ,0 
FORMAT(5X, 'at THE FINAL APPROX I M A J’l ON F'jLun.'J* //(oE2n.9j/j 
GOTO 8C0 

CONIINUE . . 

OPENCUMir=ND,DEVICE=' DSK') 

Kr<TTE(ND,752> 

FORMATf//5X, 'MODIFIED JACOBIAh IS ■ SINGULAR. I'Rv h uIFFS»E-lr'| 
KRJ IE(ND,7525) , . . 

FORHAKSX, 'INITIAL AP PR3X1NATI Oh . ' ) 

ASTER=3, 

'T¥P£ 339 ' 

FOPMAT(//5X, 'MODIFIED' JACORiAN iS SINGULAR. TH> a i ' 7 

TYPE 401 

FUFMAK5X, 'INITIAL AP PP jX 1 M ATiO.v . ' ) 

MIiXIT = Ml + l 
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CLOSF(UNlT=MO,Dfc:Vi:E,= '[)SK' ) 

RETURN 

END 

subroutine BACK(K'«IM,W,X,ISUU,Cnr;,l,njKUPJ 

Tnld suoroutlae ba:c-SDives tf.e first KmT:'* oi n LriannaVa^- 

rized lin.ear syster far In proved x values in terms ofniev-io.ib 
ones. 


10 

ICO 


DiMENSTUM X ( 30 ) , C 3E ( 3 , 3 1 ) 
dimension ISUIH3D) ,LJ0NUpf30,3i) 

DO 200 KK = 1 , NMIN 
K.v r KMiN-KK + 2 
KNRX = ISUB(K^-l) 

XtKMAX) = O.OE+00 
DO 100 J=KV,N 

JSnBsuOONUFtKM . J) 

X ( K AX) = X I K MAX) + rOE ( K v,-l , J SUB J * X f J Sli R j 
CUNTI:mIJE 

XtKMAX) = X(KMAX)4:jEtKM-l + 

CON’rlwUE 

PETURiv 

END 

subroutine A(lXFCNCX,y,Kj 
dimension xt30) 

COMMUfJ FO , AK , AKA, AXT , AKT.A, AKl , aKS,SXS/SC\, S'-.A , V , vG.u 

p,cxo,cso,cco,czo,:ro 

cx=x(n 

CS=X(2) 

CC=X(3) 

PC=X(4) 

CH0=CZ0-:50 

^ |_J — r* 2 ^ Q ^ g 

CriS=AKl*CC»:S/(AKA* tczo-rs) J 
AKU = HmII/( 1 + AKS/CHS + :.HS/AKI ) 

CuNri.NUE 


^ ^ n R ,i 


GO ro (1.2,3,4),K 

y=FO* (CXO-OX) /V+AMUtOX-AKitCTU 


RETURN 

y=FOTtCSO-CS)/V^A nJtOx/SXS 
RETURN 

y=FO* CCC 0 -CC)/V + AMU* 2 XtSBX + FO*Cc:Hu-CH)/v + 

AKLA* CAKtpO-CC) 

RETURN 

y=-p*SP^N*AKLA*(AK»P:-;::)/VG-PC^(vSP*V*SMX*A'H'^rA- 

SR^V»AKl.A* CAK*PC-c:) )/yo 

RETURN 

FND 


appendix a , 2 


THIS program computes 20 REfiATi 
state and^at a particular flow* rate. 


PROGRAM SOLVE; 
const 


tyre 


var 


m=4;N = 3;K1=0,1E-5;Ka = 0.31 ? 5-1 o'^6U!Jx=f7 * u; 

KS=6.33SE-4;Ri = l,tjE-5; SP- 25 . 5 ; V-iL.O,Sflx q .o, 

MATHix=array Ll..Nrl.-Nj of real; 

FMATix^array of reai> 

DCX0,0CS0,DCCt),CS,CCr PC, CX,CHS,MO:real; 

I,J;integer; ^ ^ ^ 

OBSVARrarray t? . • ^ f 1 • i . 

IHAlRlX:ar ray f l . .S » 1 • « i'* 3 real# 

(MATIMV inverts A f^ATRjX h,' OF URDEK L A^L) THE lN'VFHJ'--0 
matrix is CAr.LED C .3 

procedure maTi ?J V C A: MA IRI X; M ; 1 nt eaer ; var crMATPixl; 

var 

T #J#K;iateqer; 
xjreai; 
begin 

for I: = 1 to N do ^ 

for J:=i to N 3 o ^ _ 71.-1 
if i = J then - LJ 
else. C C 1 , J3 ; = 
for I ; =1 to N do 
begin 

X;=Atlfll; , ^ 
for J:=i to N do 

JJ :=ALr, Jj/x; 

CLl»J35 = -t3:#Ui'A 

end; 

for J;=l to N do 

""^'5:=Arj#l3; 

it to w do 

LJ,FJ ^ 

C £ j » E j : =C I J , X J - - * r j * X . 

end 

end 

end 

end; 

^ ..n T'Di ri."*; Fprivi i ii-' s viAiRiy 

tFOUR-TRl FORMS ALL POgSIbLE 3 + 3 MATRltEv^^ f^TK J ' 

'IMATRIX*. all of TiiESF. arc, STDRp It VEB J EL HA i F I C E 
EACH OF THESE MATRICES IS EACH IN v f- KT c,U v aTKl x 

tugather form 'GIMAT 3 -^original MATRI x.i 

IS KULTiPLlEn THE CORRtSPOMD 1 N o oka 

PROCEDURE FOUR-TRltTMArRiX;FrtATlX); 
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var 


bpgin 


I , J , Ii , K : 1 nt ecjer ; 
FRjWtbooiean; 

GMATHTx»GIMAr:arr3/ ti 
T^MP/rSf-iPi : matrix; 


for l:=l to M do 
begin 

FrOia/ : =f al se ; 
for J;=l to N do 
begin 

Hi • ~ -J ? 

lt(L. = I) toeo FRd>^;=true 
If 

(FR0w=tias:) tnen l:=L+] 
f or K : = i to do 


IV 1 of 


GMATRIX [I , J,KJ : = I«MRiX ir. 

end 

end; 

for I ;=l to M ao 
begin 

for J;=l to N do 

for K;=l to W do 

TE«PtJ,K.3 : = GMATrTx[ 1,J,K1 
for J;=l to M do 
begin 

for K:=l to N do 

end; 

MATlM\/(TFMP,N,rc;Mp) ; 
for J:=l to N do 
begin 

for k: = i to A’ do 
begin 

GIMATCI, J,KJ :=Tb:f'Rrj,K3 ; 

end; 

end 


r e i ; 





end; 

for I;=l to M do 
begin 

WRirELN;WRiXELMt ' ,l);WRliF.LN; 

forJ;=itoMdo 
begin 

for K;=l to M do 
begin 


end; 


end 


end 


end; 
WkITEuN 


TE'^Pl [ J,K] :=GMAXRiXLT , J,Ki -iA * . 

WHI TEt TEMPI [J,Kj ) , 


•il ; 


(GIVEN FOUR integers i,j,K ANd 1 ^'SF:LECT' FORMS A , 

“ i-.- 'Tf.iiMRjX'.Hr-HE ^,ARl>r 


(GIVEN FOUR integers i,j,K AN^ 1 , 'SELECT' 
GOT FROM THE iTrii , jTH r XTri/ AMD ITrt RD/iS OF '' 
IS CALLED WITH 'MAT' AS THE PARAMETER.] 


•i ^ j 


i '^1 


procedure select ( i o #ic,i; integer ) ; 

HATfFMATlXr 

M,p, 0 ;INTegeR; 

begin 
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6 rAi=ito 6 do 

It tM=i) or (M=J) or (M=K) or (M=L) 
t fieri 

begin 

for P:=l to M do 

MATCQ,P] ;=I MATRIX CM, P] ; 
q:=C)+1 
end; 

WRITKLNC 'OUTPUT FOR- ' , I ; ? , J ; 2 , K : ? , I. ; 2 ) ; 

WRITELM; 

fuur_tpi(mat) 

end; 

[IMATRIxri , J] GiVt:S TtiK RATc; OF CrJftUGi:: UF Itb 0 8 V tu .A^ThPuP. 
with kESPECI ro the Iru con TROOuEU i/AKrAF?bF. 1 


begin 

READ(DCXO,OcSO,DCCO) ; 
for I ;=2 to 7 do 
begin 

for J;=i to 4 Od rF AO C UrS V Ak r 1 , J 1 j ; 
CX:=0bSVAR[ 1 , 1 ] ; 
cs:=ObSVAR[i ,2] ; 

CC: = 'ltiSvAR[i,3) ; 

PC: = OcivSVAR[ j ; 

CtiS:=Kl*CC*CS/( KA+ (C7.J-CSJ ) ; 

MU: = RMlI/(l + KS/CHS*-?H5/Ki) ; 
nciSVARCl ,5] : = S?»7¥5^.X*'»0*CX; 

OoSv Arc X , bl :=-I.JGt; K 1 t:C/( c/u-CS; ) ; 
READ UN 

end; 

for 1 ;=i to 6 do 


begin 

JHArRlXLT,lj 
TMAlRlXlT, 2 J 
IMATRIX [I , 3 j 

end; 

for I;=i to 6 do 


f JBSV/^R (2 , T j -UBS VAR I 3 , T J ) / DC X J r 
( 08 5 V AR [ , T j -u BS V A P i B , T J f / t'C S j ; 

(OBSVaR[5,Tj-OBS'^ARL7,IJ )/or.cj; 


begin 

for J;=l to N do 

WRI l'E( IMAl'RiX IT , , 7 } ) 
WRIIEDN 

end; 


WRIIELF; WKITELU; 
SLLECr 11 , 2 , 3 , 4 ); 
SE7.Ep 

SELECT C 3, 4,5,6); 
SELECT ( 1,3, 5,6); 
SELECT Cl, 4,5,6); 
SELECT (2, 3,5,6); 
SEI.ECI (2, 4, 5,6); 



end 
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appendix B .3 

♦ ♦♦♦♦♦♦ t 


dimension OBSVAPCS J £ 0BSVAU6},rt(t>, JOJ ,jVivRtt>.l UO j, v''nS^ K ■) 
THE PARAMcTERS USED ARE 

FO = ELUW RATE,L/DAX.- 

AK = HFNRX's lAa r3NSrA,MT,MJh/L Al'M. 

A^R = ACID DTSSOCIArlON TONSTAMT , ROiy i,. 

AM=lN(Ii3lTiDN COliSl'A \,r, Mjl,/L. 

AKLA=0VERADb GAS TRANSFER FILM cnfclFT^/nAY. 

AKIeIQHISaTIOn CONSTAAiT OF BT C Auf’.uM ATE , *’.0 u/ b • 

AkSsSATURaTION COLSi’A VT, M JL/f,. 

SXS = nRGANl.SM XIELU. 

SCX=C02 YiELP 
SMX=Cri4 YiEbD, 

V=D1G£STER blyUlD VJLJME.Li. 

VG=v'nLnME OF GAS SPACE, L. 

HmUsMAX. .SPECIFIC GROriTri RArft:,/bAi. 
ak.t=tuxici influence rate, /uai. 

SP = COrvVERST JN F AC 1 0 R , L/ OOb . 

P=TOTAL PRESSURE, ATvi. 

CXO,rsO,CC(),CZO,ClD=TN[.'ET CO'JCEbTiATT jf S J^' iLvGAriS.,c,Ti''i"i- 

substrate, nivSsnuVEO cu?,')et CbTiOt. aj-'j ioxt: :'[? i '■r.'--. 

i)AlA F0,AR,ARA,AKI ,ARriA,A^i cAf" S,syb,SCA,S''Af r ' r , 

CxO,CSO,^CO,CZO,Cr[)/l .■0,0,024b9i,U.UlDF-4 ,1 .Uk^ so C-^. 

0 . SBiF-l , 0.02 , 47 . 0 , 4 ? .0 , 1 0 .0 ,7.0 ,c .4 , ? . 0 , I- , b . y T-»:>i , w . j. , 
0.09,0. 009, 0.1 ,0.0/ 

common FO , AK , AK a f AK I , AKL a , AKI , f j; S , SXS , sex , S ‘U . V ' V o , -1 ' . , b. ■< 1 , 

SP,P,CXO,CSO.CCO/C7-0, :.T0,CHs,aMu 

EXTEPNAL fcn 

N=6 

lx=30 

READ(5,»)TnL 

TFAlLsO 

TJMe1=0,0 

RbAU(2l»^) TL.DIIME 
REAi>(21,*j(0BSVAKiri , 1 = 1 ,b) 
on 5 1 = 1 , b 

nbvSVAl(l)=nBSVARf 1) 

M 1=10/0X1 ME 

Here TL=riME i,ag, on NE=TNCKEMF;i»T. 

TYPE ♦,TL,D1IME,M1 
CS=3PSVAU2) 

CC=0BSVAJC3) 

rZ=3PSVAl(5j 

CHS = AM»CC»CS/ C aka* (Cz-csn 

AKU=HMU/Cl+AKS/CHS+CHS/AKi) 

integration is D3nE FR3Y. O.U TO TL 

VAI.UES of CHS niiTAlNEO ARE STUPED Jw vCtiStMj. 

DO 10 T=1 ,Ml 


rS=DRSVAl 12) 
CC=3BSVA1 (3) 


CZ=3BSVA1C5) 

CHS=Ar1*CC*CS/{AKA*1CZ-CS) ) 

T T '1^ F* ** T T M E 1 4' D T 1 H fcj 

nag ROUTINE d 62EAF I NTEGRATES A S'JJKF SISTER up D i FF FRF . T 1 r- 



Ba 


Over a ramc:e, ^ith siutabi-f tntjtai 

VARlABLt: 70 RJE|<t^VARlAfM,F~SrK|J " “ ' 


EgUATlONS 
nSl.vG THE 
here D02EAF InTegHAIES 
OBSVAl IS OF SIZE 


N oElw 


'Mr* 

b r, 


1 T .‘iK 1 


r iNai • I 
FAF? '^r.Tnnu. 

A V j i T «i F ' . 


, F C .M , , i . I f A i- F i 


rAr.l. D02EAF (lIMEl , TIME, ^F,nbSVAl ,Tni, 

CAl.b OOiItTMEI ,0RSVAU 
VCHS(1)=CHS 
T1PE*,TI'1e1 

TYPE OBSVAl ( J) , J = 1 ,5) 

continue 

read C 21 ,*)FTIS1E 

M=(FT1ME-TF)/DTiME 

DO 20 Isl ,M 

CHS^VC^S^i) 

AKU=nMU/( 1 +AKS/CHs+CHS/ AKI ) 

IF CMjD(I,10) .NE. 0) GuTj faO 
CAF.L OUrCTlMFl OBSVAI ) 

TlME=TIMEl+DTiME 

Integration is done between tihei ahh tif'e,a1ih the vai-ue 
OF CHS LAGGING. CHS IS CriHSE^ AS VCI’SCl) A xS REAuIuSiTl;. 

CALL D02EAF(TiMEl , ri'^E, N OBSVhI , TuL , FC n , W , T ,n , IF aTu) 

DO 30 J=l,(Ml-U 
VCHSC J) = VCHS( j + n 
CS=DBSVAl(2) 

CC = OBSVAIO) 

CZ=0PSVA1(5) 

VCHS(MU=AK1»CC»CS/C A(CA* (CZ-CS) ) 

continue 

STOP 

END 

FCN COMPUTES THE. FUNCTION F FUP ThE VAF uFS Ijf V J. A -4F,!:.S 

contained in ORSVAR.' 
subroutine FCnCT.OBSVAR.F) 
dimension OBSVAR(b) ,F(65 

Common eo,ar,aka, aki, a.kla, aki ,aks,sxs,scx,Sv;a. v.Vo, a*' i ^ 

SP,Pi.CXO,CSD,:CO,CZO,:TOOHS,AMj 
IF (T .GT, 1,0) FOri.-O 
CX=OBSVARll) 

CS=0BSVAR(2) 

CC=0eSVARC3) 

PC=0PSVAR(4j 

CZ=DPSVARI5) 

CT=DBSVARC6) 

C110=CZ0-CS0 

(3 H s: 7, •* C 3 

F(n=F0*(CX0-CX)/V^ AMJ*CX-AKT*Ct 
F( 2)=F0*(CS0-CS)/V-AMa»CX/SXS 

F(4) = -P't^Sp*V*AKLA=FlAK*^PC-CC)/VG-F»Cl'tSP*V*SMA*AMu*LX- 
SP*V»AKLA»r AK + PC-CC) ) /.VG 
Ft5)=F0*(CZC-CZ)/V 
F(6)=F0» (CTU-CT) / V 

F(3)=FO*(CCO-CC)/V+A.MJ»CX + SCX + FO'*(CHO-Crl)/,V + Fi?i"r ( :») + 
AKr-A*(AK*PC-CC) . ' 

WRITEC24,»),(F(I),I = J ,6) 

RETURN 

END 

output subroutine, 
subroutine OUT(TIME,03SVAR) . 

CO-MMO^ ED , a|‘? AK a! aL , AKLA , AKl , AKS , SXS , SCX , S SA . v . VG, nM j , aK i , , 


% 
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5 

10 


6 

50 


,e>) 


cxo,cso,:cOrCzo,Cf:),Ciis.Mn 

kritf: U24.*UT1M5 _ 
wr?tf^( 2 2,») ,r=J f 

,, 

S^n§ = QSs’ARC4j/!i.l03 ' 

OCH 4 = SP* AR( 1) 


CU3,PH 


RETURN 

END 



